The material below is Lecture 2 (on Energy Minimization) from Drug Discovery Computing Techniques, PharmSci 175/275 at UC Irvine. Extensive materials for this course, as well as extensive background and related materials, are available on the course GitHub repository: github.com/mobleylab/drug-computing
This material is a set of slides intended for presentation with RISE as detailed in the course materials on GitHub. While it may be useful without RISE, it will also likely appear somewhat less verbose than it would if it were intended for use in written form.
Today: Energy landscapes, energy functions and energy mimimization
# Run cell if using collab
# Mount google drive to Colab Notebooks to access files
from google.colab import drive
drive.mount('/content/drive', force_remount = True)
# Run cell if using collab
# Move into directory so that files for this lecture can be accessed
%cd /content/drive/MyDrive/drug-computing/uci-pharmsci/lectures/energy_minimization/
%ls
from IPython.display import Image
Image('images/funnel.png',width = 600)
Image('images/funnel.png',width = 600)
A potential, $U(\bf{r^N})$ describes the energy as a function of the coordinates of the particles; here ${\bf r}$ is the coordinates, boldface denotes it is a vector, and the superscript denotes it is coordinates of all $N$ particles in the system.
Image("images/funnel_labeled.png",width = 600)
Image("images/funnel_labeled.png",width = 600)
Image("images/landscape_1D.png",width = 1000)
Imagine we have two particles with particle 1 having coordinates $x_1 = 1$, $y_1 = 3$, and $z_1 = 2$, and particle 2 having coordinates $x_2 = -1$, $y_2 = 3$, $z_2 = -1$. That would give us an array like this: \begin{equation} {\bf r}^N= \begin{bmatrix} 1 & 3 & 2 \\ -1 & 3 & -1 \\ \end{bmatrix} \end{equation}
import numpy as np
r_N = np.array( [[1, 3, 2], [-1, 3, -1]], float)
print('y coordinate of first particle is %.2f' % r_N[0,1])
y coordinate of first particle is 3.00
We could compute the distance between particles as $d = \sqrt{(x_1-x_2)^2 + (y_1-y_2)^2 + (z_1-z_2)^2}$
You could code that up directly in Python, or you can use array operations:
d = np.sqrt( ((r_N[0,:]-r_N[1,:])**2).sum() )
print(d)
3.605551275463989
We're using the notation $ {\bf r}^N$ as shorthand to refer to the x, y, and z positions of all of the particles in a system. This is actually an array, or a matrix, of particle coordinates:
\begin{equation} {\bf r}^N= \begin{bmatrix} 1 & 3 & 2 \\ -1 & 3 & -1 \\ \end{bmatrix} \end{equation}Each row of that matrix has the x, y, and z coordinates of an atom in the system. *We will use this concept heavily in the Energy Minimization assignment.
The force is the slope (technically, gradient):
$f_{x,i} = -\frac{\partial U({\bf r^N})}{\partial x_i}$, $f_{y,i} = -\frac{\partial U({\bf r^N})}{\partial y_i}$, $f_{z,i} = -\frac{\partial U({\bf r^N})}{\partial z_i}$
As shorthand, this may be written ${\bf f}^N = -\frac{\partial U({\bf r^N})}{\partial {\bf r^N}}$ or ${\bf f}^N = -\nabla \cdot U({\bf r^N})$ where the result, ${\bf f}^N$, is an Nx3 array (matrix, if you prefer)
If energy function is pairwise additive, can evaluate via summing individual interactions -- force on atom k is \begin{equation} {\bf f}_k = \sum_{j\neq k} \frac{ {\bf r}_{kj}}{r_{kj}} \frac{\partial}{\partial r_{kj}} U(r_{kj}) \end{equation} where ${\bf r_{kj}} = {\bf r}_j - {\bf r_k}$. Note not all force calculations are necessary: ${\bf f}_{kj} = -{\bf f}_{jk}$
For an Lennard-Jones 38 cluster:
Image("images/dgraph_lj38.png",width = 700)
See also Doye, Miller, and Wales, J. Chem. Phys. 111: 8417 (1999)
Note related "Python energy landscape explorer" (PELE)
(Image source https://pele-python.github.io/pele/disconnectivity_graph.html, CC-BY license, by the Pele authors: https://github.com/pele-python/pele/graphs/contributors)
Here's a bonus disconnectivity graph:
Image("images/Fedg.png",width = 700)
(Image source https://commons.wikimedia.org/wiki/File:Fedg.png#filelinks, by Snd0, CC-BY-SA 4.0)
Consider $U(x,y) = x^2 + (y-1)^2$ for a single particle in a two dimensional potential. Finding $\nabla\cdot U = 0$ yields:
$2x = 0$ and $2(y-1)=0$ or $x=0$, $y=1$
simple enough
But in general, N dimensions means N coupled, potentially nonlinear equations. Consider \begin{equation} U = x^2 z^2 +x (y-1)^2 + xyz + 14y + z^3 \end{equation} Setting the derivatives to zero yields:
$0 = 2xz^2 + (y-1)^2 +yz$
$0 = 2x(y-1) + xz + 14$
$0 = 2x^2z + xy + 3z^2$
Volunteers?? It can be solved, but not fun.
And this is just for a single particle in a 3D potential, so we are typically forced to numerical minimizations, even when potential is analytic
Common problem: For some $f(x)$, find values of $x$ for which $f(x)=0$
Many equations can be re-cast this way. i.e., if you need to solve $g(x) = 3$, define $f(x) = g(x)-3$ and find $x$ such that $f(x)=0$
If $f(x)$ happens to be the force, this maps to energy minimization
As a consequence: Algorithms used for energy minimization typically have broader application to finding roots
Here we'll set up a simple function to represent an energy landscape in 1D, and play with energy minimizing on that landscape.
#Import pylab library we'll use
import scipy.optimize
#Get pylab ready for plotting in this notebook
%pylab inline
#Define a range of x values to look at in your plot
xlower = -5 #Start search at xlower
xupper = 5 #End search at xupper
#Define a starting guess for the location of the minimum
xstart = 0.01
#Create an array of x values for our plot, starting with xlower
#and running to xupper with step size 0.01
xvals = np.arange( xlower, xupper, 0.01)
Populating the interactive namespace from numpy and matplotlib
#Define and the function f we want to minimize
def f(x):
return 10*np.cos(x)+x**2-3.
#Store function values at those x values for our plot
fvals = f(xvals)
#Do our minimization ("line search" of sorts), store results to 'res'
res = scipy.optimize.minimize(f, xstart) # Apply canned minimization algorithm from scipy
#Make a plot of our function over the specified range
plot(xvals, fvals, 'k-') #Use a black line to show the function
plot(res.x, f(res.x), 'bo') #Add the identified minimum to the plot as a blue circle
plot(xstart, f(xstart), 'ro') #Add starting point as red circle
#Add axis labels
xlabel('x')
ylabel('f(x)')
Text(0, 0.5, 'f(x)')
Try adjusting the above to explore what happens if you alter the starting conditions or the energy landscape or both. You might try:
xstart) so it is slightly to the left or slightly to the right+x. Can you make it so one of the wells is a local minimum, perhaps by altering the coefficient of this term? Take ${\bf f}^N = -\frac{\partial U({\bf r}^N)}{\partial {\bf r}^N}$, then:
Image("images/steepest_descent.png",width = 700)
Repeat until minimum is found
Limitations:
Oscillates in narrow valleys; slow near minimum.
${\bf f}^N$ is the force on all of the atoms, as an array, where each row of the array is the force vector on that atom.
$U({\bf r}^N)$ is the potential energy, as a function of the positions of all of the atoms.
These use the same vector and array notation we introdued above for ${\bf r}^N$.
Image("images/steepest_descent.png",width = 700)
#Image("images/Banana-SteepDesc.gif",width = 700)
Image(open('images/Banana-SteepDesc.gif','rb').read())
(In this case, steepest ascents, but it's just the negative...)
#Image("images/gradient_ascent_contour.png",width = 700)
Image(url='https://upload.wikimedia.org/wikipedia/commons/d/db/Gradient_ascent_%28contour%29.png')
Image(url="https://upload.wikimedia.org/wikipedia/commons/6/68/Gradient_ascent_%28surface%29.png")
A line search is an efficient way to find a minimum along a particular direction
SciPy has lots of functions and tools for optimization, but it doesn't even implement steepest descents because of poor reliability.
However, the Nelder-Mead minimization method applies a downhill simplex method which also is less than ideal, so let's play around with that a bit. First, let's make a suitable landscape:
#Define and the function f we want to minimize
def f(arr):
return 10*np.cos(arr[0])+arr[0]**2-3.+arr[1]**2
#Define a range of x, y values to look at in your plot
# NOTE IF YOU WANT TO ADJUST THESE YOU NEED TO RE-RUN ALL THREE CELLS (shift-enter)
xlower, ylower = -5, -5 #Start search at xlower and yupper
xupper, yupper = 5, 5 #End search at xupper and yupper
#Define a starting guess for the location of the minimum
xstart, ystart = -1.0, 1.0
#Create an array of coordinates for our plot
xvals = np.arange( xlower, xupper, 0.01)
yvals = np.arange( ylower, yupper, 0.01)
# Make a grid of x and y values
xx, yy = np.meshgrid(xvals, yvals)
#Store function values at those x and y values for our plot
fvals = f(np.array([xx,yy]))
colors = np.linspace(0,1,len(xx))
#Create 9''x7'' figure
plt.figure(figsize=(9, 7))
#Plot the Energy Landscape with a colorbar to the side
plt.contourf(xvals, yvals, fvals, 10)
plt.colorbar()
plt.show()
res = scipy.optimize.minimize(f, np.array((xstart,ystart)), method='Nelder-Mead', options={'return_all':True})
plt.figure(figsize=(9, 7))
plt.contourf(xvals, yvals, fvals, 10)
plt.colorbar()
# Plot path of minimization
xvals = [ entry[0] for entry in res.allvecs ]
yvals = [ entry[1] for entry in res.allvecs ]
colors = np.linspace(0,1,len(xvals))
plt.scatter(xvals, yvals, c=colors, cmap='spring', s=5, marker = "o")
plt.plot(xvals, yvals, 'm-')
[<matplotlib.lines.Line2D at 0x7f81d0f475e0>]
Start with an initial direction ${\bf v}$ that is downhill, move in that direction until a minimum is reached.
Compute a new direction using ${\bf v}_i = {\bf f}^N_i + \gamma_i {\bf v}_{i-1}^N$ where $\gamma_i = \frac{({\bf f}^N_i-{\bf f}^N_{i-1}){\bf f}^N_i}{{\bf f}^N_{i-1} {\bf f}^N_{i-1}}$; note $\gamma_i$ is a scalar.
Note that by ${\bf f}^N {\bf f}^N$ we mean vector multiplication, not a matrix multiplication; that is: \begin{equation} {\bf f}^N {\bf f}^N = f^2_{x,1} + f^2_{y,1} + f^2_{z,1} + f^2_{x,2} + ... + f^2_{z, N} \end{equation}
(In Python, this can be coded by multiplying ${\bf f}\times {\bf f}$, where ${\bf f}$ are arrays, and taking the sum of the result; for normal matrix multiplication of arrays one would use a dot product)
$\gamma_i$ is designed so that the new direction is conjugate to the old direction so that the new step does not undo any of the work done by the old step (causing oscillation)
Ideally takes at most $M$ steps, where $M$ is the number of degrees of freedom, often $3N$, though in practice even small precision errors make it take longer than this.
(Image: Wikipedia, public domain. Oleg Alexandrov.)
#Define and the function f we want to minimize
def f(arr):
return 10*np.cos(arr[0])+arr[0]**2-3.+arr[1]**2
#Define a range of x, y values to look at in your plot
# NOTE IF YOU WANT TO ADJUST THESE YOU NEED TO RE-RUN ALL THREE CELLS (shift-enter)
xlower, ylower = -5, -5 #Start search at xlower and yupper
xupper, yupper = 5, 5 #End search at xupper and yupper
#Define a starting guess for the location of the minimum
xstart, ystart = -1.0, 1.0
#Create an array of coordinates for our plot
xvals = np.arange( xlower, xupper, 0.01)
yvals = np.arange( ylower, yupper, 0.01)
# Make a grid of x and y values
xx, yy = np.meshgrid(xvals, yvals)
#Store function values at those x and y values for our plot
fvals = f(np.array([xx,yy]))
colors = np.linspace(0,1,len(xx))
res = scipy.optimize.minimize(f, np.array((xstart,ystart)), method='CG', options={'return_all':True})
plt.figure(figsize=(9, 7))
plt.contourf(xvals, yvals, fvals, 10)
plt.colorbar()
# Plot path of minimization
xvals = [ entry[0] for entry in res.allvecs ]
yvals = [ entry[1] for entry in res.allvecs ]
colors = np.linspace(0,1,len(xvals))
plt.scatter(xvals, yvals, c=colors, cmap='spring', s=5, marker = "o")
plt.plot(xvals, yvals, 'm-')
[<matplotlib.lines.Line2D at 0x7f81b42114c0>]
#Define and the function f we want to minimize
def f(arr):
return 10*np.cos(arr[0])+arr[0]**2-3.+arr[1]**2
#Define a range of x, y values to look at in your plot
xlower, ylower = -5, -5 #Start search at xlower and yupper
xupper, yupper = 5, 5 #End search at xupper and yupper
#Define a starting guess for the location of the minimum
xstart, ystart = -0.1, 1.0
#Create an array of coordinates for our plot
xvals = np.arange( xlower, xupper, 0.01)
yvals = np.arange( ylower, yupper, 0.01)
xx, yy = np.meshgrid(xvals, yvals)
#Store function values at those x and y values for our plot
fvals = f(np.array([xx,yy]))
res = scipy.optimize.minimize(f, np.array((xstart,ystart)), method='BFGS', options={'return_all':True})
plt.figure(figsize=(9, 7))
plt.contourf(xvals, yvals, fvals, 10)
plt.colorbar()
xvals = [ entry[0] for entry in res.allvecs ]
yvals = [ entry[1] for entry in res.allvecs ]
colors = np.linspace(0,1,len(xvals))
plt.scatter(xvals, yvals, c=colors, cmap='spring', s=5, marker = "o")
plt.plot(xvals, yvals, 'm-')
[<matplotlib.lines.Line2D at 0x7f81b6fa4b80>]
One can simply do lots of minimizations from random starting points and find the global minimum sometimes, if there are not too many minima in total. See e.g. this local and global minima discussion.
As an exercise, you might try to find the global minimum of this function by picking random starting points and minimizing many times (for solutions, see that link):
def f(x, offset):
return -np.sinc(x-offset)
x = np.linspace(-20, 20, 100)
plt.plot(x, f(x, 5));